Critical collapse and the primordial black hole initial mass function 
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It has normally been assumed that primordial black holes (PBHs) always form with mass approx- 
imately equal to the mass contained within the horizon at that time. Recent work studying the 
application of critical phenomena in gravitational collapse to PBH formation has shown that in 
fact, at a fixed time, PBHs with a range of masses are formed. When calculating the PBH initial 
mass function it is usually assumed that all PBHs form at the same horizon mass. It is not clear, 
however, that it is consistent to consider the spread in the mass of PBHs formed at a single horizon 
mass, whilst neglecting the range of horizon masses at which PBHs can form. We use the excursion 
set formalism to compute the PBH initial mass function, allowing for PBH formation at a range of 
horizon masses, for two forms of the density perturbation spectrum. First we examine power-law 
spectra with n > 1, where PBHs form on small scales. We find that, in the limit where the number 
of PBHs formed is small enough to satisfy the observational constraints on their initial abundance, 
the mass function approaches that found by Niemeyer and Jedamzik under the assumption that all 
PBHs form at a single horizon mass. Second, we consider a flat perturbation spectrum with a spike 
at a scale corresponding to horizon mass ~ 0.5A/q, and compare the resulting PBH mass function 
with that of the MACHOs (MAssive Compact Halo Objects) detected by microlensing observations. 
The predicted mass spectrum appears significantly wider than the steeply-falling spectrum found 
observationally. 



PACS numbers: 98.80.-k astro-ph/990126£ 



I. INTRODUCTION 

There has been a lot of interest recently in numerical 
studies of gravitational collapse, carried out for a wide 
range of matter models, which appear to exhibit critical 
phenomena: self-similarity, universality and power-law 
scaling of the black hole mass (for a review and extensive 
references see e.g. 0). Consider a smooth one-parameter 
family of initial data described by a parameter p, such 
that for p > p c a black hole is formed and for p < p c no 
black hole is formed. Numerical simulations (originally 
carried out by Choptuik for the case of a massless scalar 
field Q and by Evans and Coleman ||] for a perfect fluid 
with equation of state p = p/3) show that the black hole 
mass, Mbh, scales as 

M BU <x(p-p c y , (1) 

for p ~ p c , where 7 is a universal scaling exponent which 
is independent, for a given matter model, of the choice 
of p and the initial shape of the density fluctuations. 

It has been pointed out by Niemeyer and Jedamzik Q 
that this work has an astrophysical application to primor- 
dial black hole (PBH) formation. For critical phenomena 



to occur, the initial distribution of fluctuations must be 
such that most of the collapsing fluctuations have magni- 
tude not too far above the critical magnitude for collapse. 
This situation arises in the formation of PBHs from den- 
sity perturbations in the early universe, since the distri- 
bution of large PBH forming fluctuations falls off rapidly, 
roughly as P(S) oc exp(— 8 2 ) where 6 = Sp/p. Generally 
this is not the case during astrophysical black hole for- 
mation, such as stellar collapse where mass scales such 
as the Chandrasekhar mass are introduced due to degen- 
eracy pressure and other effects violating the scale-free 
behaviour .0 

Niemeyer and Jedamzik [El have investigated the evo- 
lution of spherically-symmetric perturbations of the en- 
ergy density in unperturbed Hubble flow during radia- 



* Such physical effects may also occur on small scales in the 
case of PBH formation so that arbitrarily small mass PBHs 
will not be formed, but this does not affect the general validity 
ofEq. (§). 
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tion domination^ for three different perturbation shapes 
at horizon crossing: a gaussian overdensity, which tends 
to the background Friedmann-Robertson-Walker metric 
at infinity, a Mexican Hat function, and an unspecified 
fourth-order polynomial. In the latter two cases the in- 
ner overdensity is compensated by a surrounding under- 
dense region. The initial conditions are tuned to super- 
criticality by adjusting the amplitude of the perturba- 
tions. Their simulations confirm the existence of the mass 
scaling relation in the case of PBH formation: 



M BH = kM H (6 - tf c )" 



(2) 



where Mh is the horizon mass at the time the fluctuation 
enters the horizon, 8 is defined as the additional mass 
inside the horizon radius in units of the horizon mass 
and 7 w 0.37 for all three shapes of perturbation. For 
gaussian shaped fluctuations k = 11.9 and S c — 0.70, for 
the Mexican Hat fluctuations k = 2.85 and S c = 0.67, and 
for the fourth-order polynomial k = 2.39 and S c = 0.71, 
suggesting that 8 C ~ 0.7 for all perturbation shapes. This 
value of 8 C is about a factor of 2 larger than the value 
found analytically || by requiring that the over-dense 
region exceeds its Jeans length at the time that it stops 
expanding. 

The formation, at a fixed horizon mass, of PBHs with a 
range of masses is a significant change from the standard 
picture where it is assumed that all PBHs are formed 
with Mbh ~ Mh Niemeyer and Jedamzik have 

determined the PBH initial mass function under the as- 
sumption that all PBHs form at the same horizon mass. 
This assumption is often used; for example, with a power- 
law spectrum with spectral index greater than one it is 
assumed that the vast majority of the black hole for- 
mation occurs at the shortest possible scale. The usual 
justification is that PBH formation is extremely sensitive 
to the amplitude of the perturbations, and hence inher- 
its a strong dependence on scale even if the variation of 
the spectrum is weak. It is not clear, however, that it 
is consistent to consider the spread in the mass of PBHs 
formed at a single horizon mass, whilst neglecting the 
range of horizon masses at which the PBHs are formed. 
Investigation of this issue is the main purpose of this pa- 
per. In the end, we shall in fact find that this assumption 
works very well in the cases of astrophysical interest. 

Hcuristically, to find the mass of a PBH formed at 
a given point in space we need to smooth the density 
field at that point, <5(x), on a range of mass scales, to 
produce 8(M), where M is the smoothing scale, and then 
evaluate each 8(M) at the time when that scale crosses 
the horizon, M = M H . If 8{M = M H ) > S c then a PBH is 



^Realistic PBH formation differs from the study of perfect 
fluid collapse in Ref. Q , where the initial data was embedded 
in an asymptotically flat space-time, in that the background 
spacetime is Friedmann-Robertson-Walker. 



formed, with mass given by Eq. (0), at that horizon mass. 
However, that condition will be satisfied for a range of 
smoothing masses; the actual mass of the PBH formed at 
the point x is given by the largest value of Mbh found as 
the smoothing scale is decreased, which will not be the 
largest smoothing scale giving a density contrast above 
the threshold due to the dependence of the mass on (8 — 
8 C ). To find the optimal smoothing mass, we use the 
excursion set formalism || which was introduced in large- 
scale structure studies to determine the mass function, 
merger rates and clustering of collapsed objects e.g. [||,[|. 



II. THE EXCURSION SET FORMALISM 

A. Standard formalism 

To examine the density perturbations on a given scale 
we must smooth the density field, as described above, 
using a window function W(Rf , |x — x'|) with radius Rf . 
The density contrast is defined as 8{x) = (p(x) — p)/p, 
and the smoothed version is given by: 

/>oo 

8(R { , x) = / W(R { , |x - x'|)<5(x')d 3 x' 
Jo 

1 f°° 

— — — 7^- / lF(fc J R f )J(k)exp(-ik.x)d 3 k, (3) 
( 27r ) Jo 

where W(kR{) and 5(k) are the Fourier transform of the 
window function and the unsmoothed density contrast 
respectively, with k = |k|. The variance of the mass 
distribution a 2 (Rf), as defined in p0|, is given by 



a 2 (Rf) = J W{kRi)V 6 (k)^ 



where T > s(k) is the spectrum of 8 



?,(fc) = -(|i(k)| 2 ) 



(4) 



(5) 



The effect of varying Rf, at fixed time, can be found 
by differentiating Eq. (§) M: 



d8{Rf,x) 
dRt 



1 



(2tt) 3 



x exp (— ik.x)dk = rj(R{ , x) 



This has the form of a Langevin equation; the change 
in (5(i?f,x) when Rf is changed is given in terms of a 
stochastic force ?y(i?f,x), which depends on the form of 
the window function used. 

A particularly good choice of window function if one 
wants to make analytical progress is the sharp fc-space 
window function 



W(kRf) = 6(1 - kRf) . 



(7) 
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Its strength is that the only effect of decreasing the 
smoothing radius is to add new Fourier modes of the 
unsmoothed density contrast <5(x) to the integral for 
a 2 . For gaussian perturbations, as we will be assum- 
ing throughout |] these new modes are uncorrelated with 
the ones already included in the integral on a larger 
smoothing scale. The change in 8(Rf,x) caused by the 
new modes is therefore independent of its value on the 
larger smoothing scale, leading to a random walk without 
memory. For other choices of window function, includ- 
ing the commonly-used top-hat and gaussian forms, this 
nice property doesn't hold since changing Rf alters the 
contribution to tr 2 from modes of all wavenumbers. That 
said, although the sharp fc-space window function is cal- 
culationally advantageous, one wouldn't expect physical 
results to be all that dependent on the choice of smooth- 
ing. 

With the sharp fc-space window function, Eq. (||) sim- 
plifies to 



d6(x, A) 
<9A 



C(A), 



(8) 



independent of position, where A = a 2 (Rf) is a pseudo- 
time variable and 

(C(A 1 )C(A 2 ))=<5 D (A 1 -A 2 ), (9) 

the right-hand side being the Dirac delta function. 
Eq. (|8|) can be integrated to give 

<J(A + 7) - 5(A) - V7G (10) 

where G is a gaussian random variable with mean zero 
and unit variance, and since this equation is exact the 
step size, 7, can be chosen freely. Stochastic processes 
which are governed by this equation, such as self diffu- 
sion in a hard sphere gas, are known as Wiener processes. 
We can think of the values of 8(A) produced as A is in- 
creased (or equivalently R{ decreased) as mapping out a 
trajectory, analogous to the path of a self-diffusing par- 
ticle. Each trajectory has the initial condition 8(0) = 
since in the limit that Rf — > 00, <5(i?f,x) — > by defini- 
tion of the mean density. 

Using Eq. ( |Io| ) and a random number generator, we can 
generate an ensemble of trajectories each representing the 
variation, with smoothing scale, of the smoothed density 
field at a different point in space. The trajectories of 8(A) 
produced for a chosen range of A values are independent 
of the form of the power spectrum. To relate each value 
of A to a mass scale, M, we need to choose a form for 
the power spectrum and then use Eq. (Q) to find the 
relationship between A and M. 



"'"The assumption of gaussianity has been challenged for large 
PBH-forming fluctuations p"l| ; however it is reasonable to 
maintain this assumption in order to assess the effect of crit- 
ical collapse on the PBH mass function. 



B. Application to PBH formation 

The trajectories give the variation of 8(M) at fixed 
time; however, the condition for PBH formation on a 
given scale, S(M) > S C: and Eq. (|^) which gives the mass 
of the PBH produced, apply when that scale crosses the 
horizon, i.e. M = Mh- We must therefore choose a 
fixed time at which to generate our trajectories, and then 
evolve each value of 8(M) forwards, or backwards, in time 
to the epoch at which that scale crosses the horizon. Dur- 
ing radiation domination perturbations which are outside 
the horizon grow as ]10[ 



8 oct (x Mh 



and therefore 



4c (M H ) = 5ft (M H ) 



Mh,d 



(11) 



(12) 



where 'he' and 'ft' denote quantities evaluated at horizon 
crossing and the chosen fixed time respectively. For PBH 
formation we require Sft(Mn) > <5 c ,ft(MH) where 



8 c ,it( M n) = 8 C 



M H 



Mi 



H.min 



(13) 



III. POWER-LAW POWER SPECTRA 

We first examine the case where the primordial power 
spectrum is a power-law 



T(S) oc k 



71 + 3 



(14) 



[equivalently P(k) oc k n , where P(k) — (|<5k| 2 ) and n is 
known as the spectral index] . Inserting this form for the 
power spectrum into Eq. (||) gives <J 2 (R) oc R~( n+3 h 

During radiation domination the mass in a comoving 
region evolves, reducing as 1/a. Ultimately we are in- 
terested in the mass associated with a given scale when 
that scale crosses the horizon; horizon crossing is given 
by R oc 1/aH oc i 1 / 2 oc M H /2 @] so that, still at fixed 
time, 



tx 2 (M H ) = (7 2 (M H , m i„) 



/ M H 



-(n+3)/2 



(15) 



where we have chosen our fixed time as the epoch imme- 
diately after the reheating period at the end of inflation 
when the horizon mass has its minimum value, MH, m in, 
which is determined by the reheat temperature, Trh, via 



M 



H.O 



RH 



To 

Tr 



cq 



3/2 



(16) 



We can obtain er(MH.min) from the mass variance on 
the present horizon scale, using the variation with M 



3 



of the mass variance at horizon crossing (for details see 
Ref. @): 



c(M H , m m) = o-(M H , ) 



M 



(17) 



where '0' and 'eq' denote quantities evaluated at 
the present day and matter-radiation equality, and 



= 9.5 x 10- 



using the COBE normaliza- 



tion) 
tion p2|. 

There are a number of well-known constraints on the 
abundance of PBHs over a wide range of mass scales; 
for an up-to-date review see Ref. Typically only of 
order 10~ 20 of the energy density of the universe, at the 
time that they form, can go into PBHs. This leads to the 
constraint that n must be less than about 1.25 pl|, with 
some dependence on the value of M Hjm i n . It is obviously 
not feasible to run simulations where only one trajectory 
in 10 20 produces a PBH. Using larger values of n, around 
1.3, a manageable number of trajectories will form PBHs 
but the values of &h c (Mn), and hence Mbh produced will 
typically be larger than when n ~ 1.25. We can however 
examine how the distribution of PBH masses behaves as 
n is decreased and compare its limiting behaviour with 
the mass function found analytically by Niemeyer and 
Jedamzik under the assumption that all PBHs form at 
the same horizon mass [01 



dN 



1 



l/ 7 



d (In .Mbh) 



(MmA 
(S c + (MbhAMh, 



x exp 



2a 2 



■ (18) 



where dN is the number of PBHs per logarithmic mass 
interval d(lnMBn), and a = ct(Mh). The Niemeyer and 
Jedamzik mass function (NJMF) will be valid if the dis- 
tribution of the horizon masses at which the PBHs are 
formed is close to a delta function. 

As an example we choose Mh.buh = 10 10 g, which corre- 
sponds to Trh ~3x 10 11 GeV, and in Eq. (||) we use the 
values of k and S c from Ref. p] for Mexican hat shaped 
fluctuations. Our general conclusions will however be in- 
dependent of these choices. Since the probability of PBH 
formation falls off rapidly with increasing horizon mass 
we generate trajectories using 1000 equally spaced steps 
in A, giving adequate coverage of the range of horizon 
masses at which PBHs form. Fig. [I] shows a typical PBH 
forming trajectory with n = 1.3. 

We ran simulations producing 1000 PBHs for each 
of n=1.310, 1.305, 1.300, and 1.295 and in each case 
found the distributions (smoothed with a gaussian to re- 
move the effects of discreteness from the finite number 
of PBHs) of the PBHs masses and the horizon masses 
at which they formed. We denote these as x(Mbh) 
and x(Mh) respectively, so that the number of PBHs 
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FIG. 1. A typical PBH forming trajectory resulting from 
a power-law spectrum with n = 1.3. The dotted points are 
the values of <5(Mh), the long-dashed line shows the thresh- 
old, at the fixed time for PBH formation 5 c ,ft(Afa), and the 
short-dashed line shows the typical size of perturbations on 
each scale, ct(Mh)- 



per logarithmic mass interval produced in our simula- 
tions is given by x(^BH)d(ln M). The upper panel of 
Fig. |^ shows the smoothed mass distributions x(Mbh)j 
along with the Niemeyer and Jedamzik mass function 
(smoothed on the same scale to allow direct comparison) 
evaluated for Mh = 10 10 g and a = 0.032, the value of 
cr(10 10 g) for n = 1.23 (the maximum value of n allowed 
for MH,min = 10 10 g by the observational constraints). 
We normalize all mass distributions to 1 at their peak to 
facilitate comparison. As n is decreased, the PBH mass 
function tends towards that found by assuming that all 
PBHs form at a single horizon mass. We can test this 
further by examining the horizon masses at which the 
PBHs form, rather than their actual masses. The lower 
panel shows the smoothed distributions of the horizon 
masses at which the PBHs formed, plotted along with 
a delta function centered on Mh = 10 10 g and smoothed 
on the same scale. As n is decreased the horizon mass 
distribution narrows. 

We therefore conclude that for power-law spectra, with 
slope satisfying the constraints from the observational 
limits on PBH abundance, the assumption that all PBHs 
form at the same horizon mass is a good approximation. 
Yokoyama ju[ has re-evaluated the constraints on the 
initial mass fraction of PBHs, taking account of the mass 
scaling relation under this assumption. He found that 
the constraints on the initial abundance of PBHs are un- 
changed if formation occurs at Mh < 5 x 10 14 g. The 
limit on the abundance of PBHs with Mbh = 5 x 10 14 g, 
which are evaporating today, coming from the 7-ray 
background, now leads to stronger limits on the initial 
abundance if 5 x 10 14 g < Mh < 10 17 g, as for horizon 
masses in this range the low-mass tail of the distribution 
has a significant population light enough to be evaporat- 
ing today. 
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FIG. 2. The smoothed distributions of the actual PBH 
masses, x(Mbh)> (top panel) and the horizon masses at their 
formation, ^(Mh), (lower panel). From top to bottom on 
the right-hand side of each diagram, the solid lines are for 
n = 1.310, 1.305, 1.300 and 1.295. In the top panel the dot- 
ted line shows the Niemeyer and Jedamzik mass function, 
evaluated for <r(10 10 g) = 0.032 and smoothed on the same 
scale, while in the bottom panel the dotted line is a smoothed 
delta function, centered on Mh = 10 10 g. 

IV. SPIKY POWER SPECTRA 

Observations of microlcnsing of stars in the Large Mag- 
ellanic Cloud |L5| appear to show that a large fraction of 
the halo of the galaxy is in the form of MAssive Compact 
Halo Objects (MACHOs) with masses of M ~ O.5M . It 
has been proposed that MACHOs may be PBHs formed 
either during the QCD epoch (ll|[l7j because of the re- 
duction in pressure forces at that time, or due to a spike 
in the primordial density perturbation spectrum Q at 
the scale corresponding to Mh ~ 10 33 g. 

The mass distribution of the MACHOs is steeply 
peaked. In Ref. |L5| the observed microlensing events 
are fitted with power-law mass functions: 



V-(M) = AM a (M min < M < M max ) 



(19) 



lo£ 





o( M h/ 



0.5 



10 33 g) 



FIG. 3. The form of <j(Mh) produced by a spike in the 



density perturbation spectrum with C 
E = 2 x 10 9 Mpc~\ 



1 x 10 10 Mpc _1 and 



= (otherwise) , 

where ip(M)dM is the mass fraction of MACHOs be- 
tween M and M + dM and A is determined by the mass 
fraction of the halo in the form of MACHOs. The maxi- 
mum likelihood fit (for their 8 event sample) is found for 
a = -3.9 and M mln = 0.30M Q |] 

The spread in the masses of PBHs formed at a sin- 
gle horizon mass, found when critical collapse is taken 
into account, raises the question as to whether or not 
it is possible to produce a population of PBHs with a 
sufficiently sharply peaked mass distribution. There are 
several models of inflation which produce a spike in the 
power spectrum on a particular scale jl8| . For generality, 
we will take the spectrum of the curvature perturbations, 
Vk, which is defined, analogously to Vs, as [Eo| 



Vn= 2^ (l ^ (k)l> 



(20) 



to have the form of a gaussian spike at a scale fc c = 
fce q (M H ,cq/MH,c) - 5 = 4.25 x 10 10 Mpc"\ corresponding 
to Mh, c = 10 33 g, with variable amplitude and width, 
E, superimposed on a flat spectrum normalized to the 
COBE data on the present horizon scale. The COBE 
normalization gives Viz = 2.28 x 1CP 9 [[l0| so our spec- 
trum is 



Vn = 2.28 x 10~ 9 



( (k-k c ) 2 




(21) 



In practice the results from taking a spectrum of this 
form are quite general, as the bulk of the PBH production 



§ For a subset of 6 events (chosen to exclude a binary lensing 
event whilst maintaining the mean duration) the most likely 
mass distribution is a delta function. 
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FIG. 4. A PBH forming trajectory produced when Vn has 
a gaussian spike with E = 10 9 Mpc _1 and C = 10 10 Mpc _1 . 
The amplitude jumps dramatically when the spike is en- 
countered. The dotted points are the values of <5(Mh), the 
long-dashed line shows the threshold 6 c ,ft for PBH formation, 
and the short-dashed line shows the typical size of perturba- 
tions on each scale ct(Mh)- 

occurs right at the peak of the power spectrum. A Taylor 
expansion in the vicinity of the peak requires only the 
amplitude and curvature, and the two parameters of the 
gaussian can reproduce an arbitrary function near the 
peak. Our results will therefore prove extremely general. 

With this more complicated spectrum, in order to ob- 
tain A(M) = a 2 (M) we have to numerically integrate 
Eq. (||) for each Vu we choose, where during radiation 
domination 

so that 

Fig. 3 shows o-(Mh) when C = 1 x 10 10 Mpc -1 and E = 
2 x K^Mpc" 1 . 

We have run simulations for a range of values of C 
and E producing, in each case, 1000 PBHs, choosing the 
fixed time at which the trajectories are produced as the 
epoch when Mh = 10 33 g, corresponding to half a solar 
mass. Density perturbation spectra with the same value 
of C but different E have the same value of <t(Mh) on 
small scales (Mh < 10 33 g), decreasing as C is increased. 
Once again it is not feasible to use spectra which pro- 
duce PBH abundances consistent with the observational 
limits. We therefore choose values of C and E which 
lead to larger values of <t(Mh) on small scales (~ 0.2) 
and examine the behaviour of the PBH mass and horizon 
mass distributions at formation as C and E are varied. 
Fig. |I] shows an example of a PBH forming trajectory 
with E = 10 9 Mpc _1 and C = 10 10 Mpc" 1 . 
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FIG. 5. The smoothed distribution of the actual masses, 
x(Mbh), (upper panel) and horizon masses, x(Mh), (lower 
panel) of PBHs formed due to a spike in the density pertur- 
bation spectrum with G = 1 X 10 10 Mpc -1 . From right to left 
on the right-hand side of each diagram, the solid lines show 
E = 1 x 10 8 , 5 x 10 8 , 1 x 10 9 , 1.5 x 10 9 and 2 x 10 9 Mpc~\ The 
dotted line in the upper panel shows the smoothed Niemeyer 
and Jedamzik mass function evaluated for <r(10 33 g) = 0.053, 
and in the lower panel a smoothed delta function centered on 
M H = 10 33 g. 

The upper panel of Fig. |5|shows the smoothed distribu- 
tion of the PBH masses, x(Mbh), for C = 1 x 10 10 Mpc -1 
and a range of values of E. For comparison we also plot 
the smoothed NJMF evaluated for tx(M H = 10 33 g) = 

0. 053, the largest value compatible with the requirement 
that the present-day density of PBHs is consistent with 
the present-day age and expansion rate of the universe, 

1. e. SIpbh.o < 1- For all the values of E we used, the 
PBH mass function has the same shape as the NJMF, but 
it is centered on larger values of Mbh, shifting towards 
smaller values of Mbh as E is increased. The lower panel 
shows the smoothed distribution of the horizon masses 
at which the PBHs form, x(-^h), for the same set of C 
and E values, along with a smoothed delta function cen- 
tered on Mh = 1 x 10 33 g. The horizon mass distributions 



G 



-0.5 0.5 

log 10 (M BH /l033g) 



0.2 



0.4 



0.6 



O.f 



A 




-0.05 0.05 

log 10 (M H /10»g) 

FIG. 6. The smoothed distributions of the actual 
masses, x(Mbh)j (upper panel) and horizon masses, x(Mh) 
(lower panel) of PBHs formed due to a spike in the den- 
sity perturbation spectrum with E = 1 x 10 9 Mpc^ 1 and 
(from right to left on the right-hand side of each dia- 
gram) C = 1.75 x 10 10 ,1.5 x 10 10 ,1.25 x 10 10 ,1 x 10 10 
and 7.5 x 10 9 Mpc _1 . The dotted line in the upper panel 
shows the Niemeyer and Jedamzik mass function evaluated for 
cr(10 33 g) = 0.053, and in the lower panel shows a smoothed 
delta function centered on Mh = 10 33 g. 



have the same shape as the delta function but are cen- 
tered at smaller Mh. As £ is decreased the centre of the 
distribution tends towards Mh = 10 33 g. 

Fig. H shows the same distributions for £ = 1 X 
10 9 Mpc~ 1 and a range of C values. As C is decreased 
the centre of the PBH mass function tends to that of 
the NJMF. The distribution of the horizon masses at 
which the PBHs form is independent of C and centered 
at M H < 10 33 g. 

The horizon masses at which the PBHs form are typi- 
cally smaller than 10 33 g, since whilst V-r has a spike cen- 
tered at this scale, ct(Mh) carries on increasing, as Mh 
is decreased, until Vn C Con that scale. By increasing 
the scale on which the spike is centered we could tune 
the horizon mass distribution to be precisely centered at 



FIG. 7. The observed distribution of the MACHOs, given 
by Eq. (|l!]), is shown as the solid line and the Niemeyer and 
Jedamzik mass function, given by Eq. ([l8|), as the dashed 
line. We fix M H = 2.05 x 10 33 g so that the mass functions 
have the same mean. The curves are chosen to have the same 
integrated mass fraction and the same mean. 



M H = 10 33 g. 

The PBH masses are typically larger than given by 
the NJMF with <j(10 33 )g = 0.053, despite being formed 
at smaller Mh- This is because our spiky density pertur- 
bation spectra all have <t(Mh) ~ 0.2, in order to produce 
a larger number of PBHs than is allowed by observations 
(as is discussed above), leading to larger values of 5(M) 
and hence Mbh- 

Since the spread of Mh at which the PBHs form will be 
extremely small if the amplitude of the spike, and hence 
<t(M), is reduced so that the number of PBHs produced 
satisfies the constraint ^(pbh.o) < 1) then the PBH mass 
function will be the same, to a good approximation, as 
that found by assuming that all the PBHs form at a sin- 
gle horizon mass. We therefore compare, in Fig. [?|, the 
observed MACHO mass function, ip(M), with that de- 
termined by Niemeyer and Jedamzik, Eq. (|l8|). Defining 
-0(M)d(lnM) to be the mass fraction in a logarithmic 
mass interval d(lnM), then ip(M) = ip(M)M. Since 
dAVd(lnM) cx ij>/M, then t/»(M B h) oc diV/d(lnM B H)- 
That is to say, the mass fraction per unit mass interval is 
just proportional to the number density per logarithmic 
mass interval. We adjust the horizon mass at which the 
PBHs are assumed to form at to Mh = 2.05 x 10 33 g, so 
that the PBH mass distribution has the same mean mass 
as the MACHO distribution. The mass distributions are 
normalized to give the same total mass fraction when 
integrated. 

The PBH mass distribution is much broader than that 
fitted to the observed microlensing events; its full width 
at half maximum is six times larger (assuming the mass 
functions are normalized to the same mean). Because 
the derivation assumes all black holes form at the same 
horizon mass, it is the narrowest possible mass func- 
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tion which can arise if the critical collapse hypothesis 
is correct. The mass function is almost symmetrical, ex- 
tending to Mbh "C Mh- Currently there is a robust 
limit, due to the absence of short-duration microlens- 
ing events, that less than 20 percent of the dark mat- 
ter in the halo of the galaxy can be in low mass MA- 
CHOs (10- 4 Af o < M < 0.03M Q ) ||, but the PBH 
mass function is easily consistent with that. Whether or 
not the width of the predicted mass function is in agree- 
ment with the observed microlensing events, which so far 
have been fitted with much narrower functions, is a much 
more pertinent question. As the duration of microlens- 
ing searches increases, and further microlensing events 
are observed, the MACHO mass function will be deter- 
mined much more accurately, and if the MACHOs are 
PBHs then the true MACHO mass function will be con- 
siderable wider than the sharply-peaked mass functions 
which have been fitted to the current data. 



V. CONCLUSIONS 

Critical gravitational collapse, and in particular the 
scaling relation for the black hole mass, has an astro- 
physical application to PBH formation. We have used 
the excursion set formalism to determine the PBH mass 
function, when formation on a range of mass scales is 
taken account of, for two different types of power spec- 
tra. The first is power-law density perturbation spectra, 
and the second is flat spectra with a spike on a given 
scale. In both cases we find that as the parameters are 
adjusted so that the abundance of PBHs decreases, the 
mass function tends towards that found by Niemeyer and 
Jedamzik under the assumption that all the PBHs form 
at a single horizon mass. There are tight observational 
limits on the abundance of PBHs, and for power spectra 
which satisfy these constraints the assumption that all 
PBHs form at the same horizon mass is a good approxi- 
mation. 

As an application, we then compared the Niemeyer 
and Jedamzik PBH mass function with that of the 
MACHOs found from microlensing observations. The 
PBH mass distribution is considerably broader than the 
sharply-peaked simple mass distributions which have so 
far been fitted to the observed microlensing events. If 
the MACHOs are indeed PBHs, then future microlens- 
ing searches should unveil the broad spread of the PBH 
mass distribution. 
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